Continuous and Discrete Adjoints to the Euler 
Equations for Fluids * 

Frederic Alauzet^ Olivier Pironneau * 
January 19, 2013 

Abstract 

Adjoints are used in optimization to speed-up computations, simplify 
optimality conditions or compute sensitivities. Because time is reversed 
in adjoint equations with first order time derivatives, boundary conditions 
and transmission conditions through shocks can be difficult to understand. 
In this article we analyze the adjoint equations that arise in the context 
of compressible flows governed by the Euler equations of fluid dynamics. 
We show that the continuous adjoints and the discrete adjoints computed 
by automatic differentiation agree numerically; in particular the adjoint is 
found to be continuous at the shocks and usually discontinuous at contact 
discontinuities by both. 

Introduction 

In optimization adjoints greatly speed-up computations; the technique has been 
used extensively in CFD for optimal control problems and shape optimization 
(see [21 [201 H21 IIZ1 HI] and their bibliographies) . Because time is reversed in 
adjoint equations and because convection terms operate in the opposite direc- 
tion, boundary conditions and transmission conditions through shocks can be 
difficult to understand. Automatic differentiation in reverse mode [10], [13] auto- 
matically generates the adjoint equations so the problem does not arise because 
boundary conditions are set by the adjoint generator. But the problem then is 
to understand and make sure that there is a limit when the mesh size tends to 
zero and that this limit agrees with the continuous solution. 

Answering these questions is important for design of supersonic airplanes 
because shocks are involved; a good case study is the design of airplaines with 
the least sonic boom at ground level. Several investigations have been done 
already (see B. Mohammadi et al [50], S. Kim et al [TS], A. Loseille et al [T§] . 
M. Nemec [21] and the bibliography therein) and automatic differentiation is 
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known to work also in this case. Another application is the optimization of 
"goal oriented meshes" [3] which minimize the numerical error to compute a 
functional like drag or lift [19]. 

Finite volume schemes for the Euler equations work also for the linearized 
Euler equations, and so should work equally well on the differentiated Euler sys- 
tem. However the variables of the differentiated Euler equations are much less 
regular than those those of the Euler equations: being derivatives of discontinu- 
ous functions the sensitivities of the flow variables have Dirac singularities and 
it is not at all clear that a finite volume method is capable of computing such 
solutions; thus the discrete adjoint systems derived by automatic differentiation 
in reverse mode are based on dubious assumptions. In other words we may 
have a convergence problem: does the discrete adjoint state converge to some 
continuous state? Is it solution of a continuous adjoint equations; and since 
there may not be a unique way to set up an adjoint for a given optimization 
problem (see Giles[7]), which adjoint is it? 

All these questions have been addressed in the case of scalar conservation laws 
in one dimension of space by Ulbrich et al. [37J 13 • A framework was also intro- 
duced in |27j and [2] which can be extended to multidimensional vector systems 
of conservation laws but the theory of course is lacking. So we propose here 
to recall the formalism and do detailed numerical tests to check if the results 
obtained in one dimension of space extends to two and three dimensions. Our 
answers are encouraging: discrete and continuous extended calculus of varia- 
tions agree; the discrete adjoints seem to converge to the continuous ones; when 
the mesh is refined the method is stable and gives better solutions. 

1 Preliminaries 

1.1 Calculus of Variations 

Consider the following problem 

min{ J(m, a) : L(u) = /(a)} (1) 

a 

where u — > L{u) is an unbounded operator, like a partial differential system, 
from a Hilbert space U to another Hilbert space V . 

When a varies, everything varies up to higher order terms according to 

5 J = J'Ju + J'Ja, with L'(u) 5u = f'Ja. (2) 

For Problem ([lj, the adjoint state is u* G V* , the solution of L*u* = J' u , i.e. 

< L* u*,v>:=< u*,L'[u)v >= J'{u)v, Vu G U (3) 

where < •, • > denotes the duality product between the spaces and their duals, 
U and U* in the first instance and V and V* in the second. 
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From these definitions we can derive another expression for 6 J because 

J' u Su =< L* u* , Su >=< u* , L' Su >=< u* , f' a da > (4) 

and finally 

SJ=<f' a *u* + J a ,Sa>. (5) 

This is a very useful result because it implies that at the solution we have: 
f' a u* + J' a = 0, the so-called optimality conditions. It also says that an iterative 
scheme like 

a m+1 =a m -\(f a *u* + J' a y\a= a m (6) 
is a steepest descent algorithm which will converges for a small enough step size 

A e R+. 

This calculus is formal but can be justified under continuity and differentia- 
bility hypotheses of all functions and operators; the difficulty is to show that 
the higher order terms are indeed small. 

2 An example with the Euler equations for in- 
compressible flows 

2.1 The Continuous Case 

Consider an optimization problem for an incompressible inviscid Euler flow. The 
velocity of the fluid is denoted u the pressure p, the domain occupied by the 
fluid 0; it could be the 2D cross section of a nozzle with an inflow boundary r~, 
an outflow boundary T + and walls T°. One seeks for a nozzle inflow velocity a 
leading to a velocity Ud in a region of space Dcfi: 

min a { J(u) = ^ f \u — Ud\ 2 : d t u + V • (u ® u) + Vp = 0, V • u = 0, 

2 Jdx(o,t) 

Ut=o — 0, u\r- = a, u ■ n|r+ = b, u ■ n\-pa = 0} (7) 

where n is the outer normal to T := dfl = T~ U T + UT°. It is assumed that the 
integral of a ■ n on T~ is equal to the integral of b on T + . 

Calculus of variations is done in H , the space of square integrable divergence 
free functions in fi. 

SJ = (u- u d )6u; / v(d t 5u + V • (Su® u + u ® Su)) = (8) 

iflx(0,T) Jnx(0,T) 

for all v G H ; a ® b is the tensor aibj and below we will use the notation 
A: B = J^ij AjBij. Let u* e be solution of 

d t u* +u-Vu* - (Vu)u* +Vp* = l D (u-u d ), V-u*=0, M*| t=T = (9) 
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with u* = on r+, and u* ■ n = on T° U . 

Notice that V(m* • it) = (Vu*)u + (Vu)u*; so by integrating Q multiplied 
by Su and performing an integration by parts : 



SJ = I Su- {d t u* + u ■ Vu* + (Vu> + V(p* - u* ■ u)) 

JQx{0,T) 

(u* ■ dtSu + V(5u : u ® u* + u* V • (Su ® u)) 

Ox(0,T) 

<5u ■ (n • u u* + n ■ u* u + (p* — u* ■ u)n) (10) 

rx(o,T) 

The middle integral reduces to the integral of —u*VSp because of ^ . Finally, 
as u* ■ n = on T, 

SJ = (a ■ n u* ■ 5a + (p* — a ■ u*)Sa ■ n). (11) 

Jr-x(o,T) 

2.2 The Discrete Case 

Let us use a finite element approximation for u of degree 2 on a triangulation 
l~h and of degree 1 for p. With the characteristic-Galerkin method [23] at each 
time step one has to solve: Vv h e Vfi, Vg/, € V^ 1 , with v/Jr- = Oj^ft ' n|r = 0, 



/ (^T+Vp m+1 H ~ J q,y-K l+1 = J -u%(x - Stu%(x))v h (x)dx 
with K = {w h e C°(fi h ) : iu h | T € P^Tfc) for all triangles T fe e 7^} (12) 



It is not hard to prove (see Appendix A) that 



5 J = St 



« ' V*a fc ) + <5a h ■ « m • V^ 1 ) - Pl m V ■ <fa A ) (13) 



with the adjoint state defined by 

Sfit + Wqh)<m ~ / prv ' Uh = J D {u ™ +1 Ud)Vh 



i 



u* h m+1 (y + Stu'^(x))v h (y)dy- I v h u* h 



* m+1 



,m+l 



(y)dz 



<(T) = 0, u* h \ r + = 0, u* h ■ n\i 



0. 



(14) 



for all u/j, with t>/i|r+ = 0, Vh • Ti|r = 0. In view of the fact that the integral 
is on the boundaries triangles only and that the gradients are differences from 
the value of the functions on the boundary and zero on the opposite side of the 
triangles we see that (74) is a discretization of (11 1. 
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Figure 1: Level lines of the adjoint equation to problem In the bottom 

part of the channel u — (2, 0) T and in the upper part u = (1, 0) T . The criteria 
involves an integral of (ui — 3) 2 in the box shown in the middle of the channel. 
The plot shows the horizontal component of the adjoint state in duality with u: 
it is discontinuous on 3 horizontal lines left of the box and in the box. The figure 
on the right is a front 3D view of the same. 



2.3 Contact Discontinuities: Vortex Sheet 

What happens when the flow is discontinuous as in the case of contact discon- 
tinuities? 

As an illustration we considered Q, = (0, 7) x (0, 2) with u(0) = [l + l„<x, 0] T , 
a(y) = u(0)\ x= q and other boundary conditions chosen so that p = on the 
boundary, with u d = (3,0) T and D = (4,5) x (0.7, f. 3). Although there is 
an obvious analytical solution (the solution is stationary : u = u(0),p = 0) the 



adjoint has been computed numerically by solving ( 1 4 1 with the P 2 /P 1 element; 
the numerical test shows that the adjoint is discontinuous (fig [I]). The scheme 
handles the situation perfectly, without oscillation. 



2.4 Conclusion 

From this example we see that both the continuous and the discrete problems 
can be derived by calculus of variations; however the computations are formal 
and although the discrete formulae seem to be a discretization of the continuous 
one it does not appear easy to prove rigorously. 

Notice also that by a deriving a discrete adjoint the upwinding term for the 
adjoint equation is naturally in the right direction. 



3 Generalization to Shocks 

In the presence of shocks, several problems arise with formidable mathematical 
difficulties: existence and uniqueness of solution, well posedness of the adjoint 
equations, continuity of these with respect to boundary data, etc. But there is 
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one difficulty which even a mindless practitioner cannot ignore: shocks! Indeed 
even the formal calculation no longer makes sense in the presence of shocks. 

For example if the density p and the velocity u are discontinuous, but with 
pu continuous, it is meaningless to write 

6 (pu) = uS p + pSu (15) 

because Sp and 6u contain Dirac singularities which cannot be multiplied by a 
discontinuous function if both singularities are at the same point. 

Indeed let H(x) be the Heaviside function (zero when x < 0, 1 otherwise); 
assume that p(x) has a shock at x = x s (a) where a is a parameter: 

p = p- + {p+ - p-)H(x-x s ) (16) 

Assuming that not only the position of the shock x s but also p ± depends on the 
parameter a then by differentiation 

8p = p~' 8a + (p +> — p~')H(x — x s )Sa — x s '(p + — p~)5(x — x s )Sa 

= p' 5a — x' s [p]S(x — x s )5a (17) 

where primes denote pointwise derivatives with respect to a; this is because the 
derivative of H is the Dirac mass 5. So the Dirac-singularity of Sp is at the same 
point as the discontinuity of u. 



3.1 Extended Calculus of Variations 



Although (151 is not allowed, yet there is a way out: 

5{pu) = (8p)u + pSu (18) 
where p (and u) stands for the mean value of p at all points: 

P( x ) = \{pi x+ ) +P(x~)) 

This idea has be used in [2H1 12] to extend calculus of variation to discontinuous 
functions. Indeed when pu is continuous across the shock, there is no Dirac 
mass on the left hand side of (18); let us show that there is no Dirac mass on 
the right either. 

Let [u] denote the jump across the shock, i.e. the value on the right side of 
the shock minus the value on the left side, then on the right hand side of ( 18 1 
the weight of the Dirac mass is 

Sa 

-[p]x' s 6au- [u]x' s 5ap = -x' s — {(p + - p~)(u + + u~) + (u + - u~){p + + p~)) 
= -x' s Sa[pu] =0 " (19) 

Even if [pu] ^ the weights on the Dirac masses are the same on both sides, so 
( 18 ) is true pointwise and singularitywise at the shock. 
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More complex functions need to be decomposed into elementary products; for 
instance to compute the extended variation of p -1 one writes 

1 = ~ = -Sp + p5(-) J(-) = -iijp (20) 
p p p p pp 

Similarly p" 1 ?/ 2 being (p _1 )u with v — uu, the mean value can be computed 
as such. For non-rational functions f(p) which cannot be decomposed into ele- 
mentary products (like ship), the overline operator on derivatives is the Volpert 
ratio: 



Definition 



f'{p) ■= 



at shocks, 



[Hp)] 
Ip] 

/'(p) at regular points 



(21) 



This is consistent with (17) 



[/] 



Sf(p) = f'(p(x))6p(x) -[f]x' s S(x- x s ) = f'(p)Sp\ X jt Xa - -rj[p]x' a 6(x - x s ) 
= f'(p)Sp\x* x ,+7(p)6p\x=x. =7{p)Sp ' (22) 

For vector or matrix valued functions one should use the identity 



[F(u)] 



F„(su+ + (1 - s)u-)ds 



(23) 



and the Volpert ratio on each component. 



4 Shocks and the Adjoint of Burgers' Equation 

4.1 Analytical check 

First consider the following example: 

Given a parameter a e R, consider the entropy solution of 

u 2 

d t u + d x ( — ) = u(x,0) = (l + a){l-H(x)) (24) 

The solution is explicit and its derivative with respect to a can be computed 
analytically 

u=(l + a)(l-H(x- 1 ± a t)), 

u' a = l-H(x- l -±^t) + ^t8(x - !±^t) (25) 
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Let 



J(a) :=ir <Tf = 1(1 + a) 2 ((l + a)T 1) 



(26) 



Then J'(0) = |(3T-2). 

To recover this by the extended calculus of variation, one proceeds as follows: 



J ' = I uu 'aW = / U*u' a \ = I U*(x,0) 

with u' solution of 



(27) 



cV a + d x {uu' a ) = 0, <(0) = 1 - H(x) 



and u* defined by 



Therefore 



d t u* + ud x u* = 0, u*\t = J' u ~ u\t^\ 



u*[x,Q) = 1-\h{x + %)-\h{x-%) 



(28) 



(29) 



and (27) and (26) agree. 



Notice that u* is integrated backward in time; it is constant on the char- 
acteristics, lines of slope 1 + a on the left of the shock and of slope on the 
right. So the characteristics do not define u* everywhere: there is a "shadow" 
triangle where no upgoing characteristics cross the line t = T . However thanks 
to the overline on u in ( 28 ) there is one characteristic of slope \ (1 + a) on which 
u* = u(T), the shock line in fact. Once u* known on the shock line then any 
point in the shadow triangle has an upgoing characteristic which crosses the 
shock and so u* is also known there, constant in this example and equal to 

Notice that u* is continuous across the shock and discontinuous at the char- 
acteristics that meet the shock at t = T . 



4.2 Numerical check 

Whether a numerical scheme will be capable of such a reconstruction is a matter 
of wonder. 

Consider another example with T = 2 for which one seeks to compute the 
sensitivity with respect to a at a = 0: 



1 



2 



,2 



J ( a ) = o / \u(x,T)\ 2 dx : d t u + d x — = 0, w| = - min(atan(a; + a), 0) (30) 



R+ 



Although there is no shock at time zero, one develops later. The following 
conservative finite difference scheme has be used 



St 2Sx K ^ ' 2Sx 
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with Si = 1 if uf + uf_ x > 0, otherwise (31) 
When J is approximated by \ J2i- X >o l^fH 2 ' the discrete Adjoint to (311 is, 



j t 1- i Si - u i+1 s l+1 + u - Si) - u i (1 - s l+1 )) = 



*M M 

u , = u, 



(32) 



with Si 



SJ 



sign(w™ + u™!). Then 



E/ *0 . °' ; 0/ *0 *0 *0 /-, \ *0/ 1 nnc Or 



i = l 



On this simple example the adjoint u* seems to be correctly calculated by 
the conservative scheme (see Figure (2J). 

The results of Table [l] show that the scheme computes J correctly. Auto- 
matic differentiation of the computer program also gives the correct ^(0). On 
figure [2] the solution and its derivative with respect to a are shown. 



J(0) 


J(8a) 


(J(Sa) - J(0))/(Sa) 


J'a(P) 


A J' a (0) 


0.195009 


0.190057 


-0.4952 


-0.492863 


-0.492863 



Table 1: Scheme (31 1 gives the correct result for J and J' a . A computation 
of J' a (0) has also been done by finite difference with Sa = 0.01 then with the 
discrete gradient (32) using the discrete adjoint and finally with Automatic 
Differentiation : A J' a (0). 




■2 -1.5 -1 -0.5 0.5 1 -2 -1.5 -1 -0.5 0.5 1 



Figure 2: Left displays at T=0.325 u (green) solution of Burgers' equation and 
the adjoint state u* (red), while u' a is shown on the right. Notice the Dirac-like 
singularity of u' a . 
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4.3 Stationary problems 

The stationary case is difficult but it can be investigated as a limit case of the 
time dependent problem. As an example consider 



J-=^l u 1 : d t u + d x (^)=0 u(x,0) = (l + a)(l -#(*)) 

flx(0,T) 



2T 

Then 



J' = i / to' with <9 t u' + d x {uv!) = u'(ar, 0) = (1 - H{x)) 

Rx{0,T) 



T 

The adjoint equation is 



d t u* + ud x u* = is, u*| T = (33) 



and then at a = 



J'= / u'(0)«*(0)= / (l-fl"(x))u*(0)= /° u*(0) 
Here too w* is continuous at the shock. 



4.4 Convergence 

For scalar conservation laws in ID the extended calculus of variation has been 
justified by S. Ulbrich [27]. He introduced the concept of shift derivatives for 
functions of bounded variations on an interval I — [a,b]. The idea is that if 
p € BV(I) it may have discontinuities at a € /, so - as explained in ( 16|17 ) 



on top of the usual variation in magnitude Sp, the variations of the position of 
the discontinuities should also be considered. The result is not a function but 
a distribution, however a Dirac mass can be approximated by a large function 



on a small support, so (17) can also be written as 



Sp = p' a 6a - x s ' a [p](H(x — x s — Sa) - H(x - x s )) (34) 

Hence Ulbrich suggests to focus on what he calls shift-variations of p which 
includes at singular points a := x s of p 

S a p := sign(5a)[p(a)]~l {a _sa- ,a+5a+) (35) 

where [p]~ (x) = — min(0, p{x + ) — p{x~)) denotes the negative part of the jump 
and Sa ± the positive and negative part of da. Thus a general shift-variation of 
p is 

dp := 5p(x) + ^25 ak p (36) 
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Then u — ► y(u) is said to be shift-differentiable at u if there exists a linear 
operator y' u and a finite set of values y' k associated with a set of singularity 
points at such that 



y(u + 6u) = y(u)+y' u (u) ■ Su + ^ y' k S ak u + o(\\5u\\) (37) 

k 



With this definition discontinuous functions can be expanded by a Taylor for- 
mula but only for a restricted class of Su. 

Furthermore when it exists the shift derivative is also the distribution deriva- 



tive such as the one computed for ( 16 ) above. 

Ulbrich[27] proves that entropy solutions of scalar conservation laws in ID 
are shift differentiable with respect to initial data under very reasonable hy- 
potheses; the framework gives a natural way to define an adjoint and the duality 
formula ^ necessary for the extended calculus of variation to hold. In an other 
article Ulbrich and Giles [9] have studied the convergence of the discrete adjoint 



to the continuous one in a simple case similar to ( 30 1 



5 Extension to the Full Euler System in 2D 

Let us apply the same formalism to the two dimensional Euler equations for 
fluids: 

d t W + V • F(W) = 0,W(0) = (38) 

where W — (p, pu, pv, pE) T is the vector of conservative variables and the tensor 
F represents the convective flux: F(W) = Fi(W) e x + F 2 (W) e y with 



F X {W) 



( pu 

pu 2 + p 
puv 

V ( P E+p)u J 



\ ( 

and F 2 (W) = 



pv 
puv 
pv 2 + p 



V ( P E + p)v ) 



Additionally we assume that W is given at time t = 0, and for bounded domains 
fl we add boundary conditions which we will discuss later at length and which 
we denote loosely by We W a d- 



We wish to measure the change S J of a functional J when something changes 
in the data (parameters) of the above problem. 

For the silent business jet problem the functional is an integral of the pressure 
on a part S of the boundary of the domain (the ground) ; for the control of the 
outflow of a scramjct J is an integral of the density on the outflow boundary. 
In both cases there exists a linear map W — > B(W) € K d and a vector b £ K d 
such that 



J = 



1 



Sx(0,T) 



\B(W) - bf 



(39) 
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The parameter (denoted by a in the previous section) could be in the shape of 
the business jet or of the scramjet but at this level it is not important. The 
extended calculus of variation on J gives : 

SJ = I WW) - b)B'SW, (40) 

JSx(0,T) 

Turning to SW we know from above that it satisfies 

d t SW + V ■ (F'{W)5W) = , 5W(0) = 0, SW - ASa e W (41) 

where Wo is the tangent space to W a d at W and ASa is the (possible) change 
of boundary condition du to the change a — » a + Sa. So let us introduce the 
adjoint equation 



d t W* +F'(W) VW* = 0, W*(T)=0. (42) 



with boundary conditionsto be chosen later. In variational form (42) is 



(d t V ■ W* - F'(W) VW*-V) + / W*(0) ■ V(0) = VV. (43) 



Let us choose V = SW, integrate by parts and use (41) 

= / {d t SW ■ W* - F'iwf^W* ■ SW) 
Jnx(o.T) 



Hence 



W* ■ (d t SW + V ■ (F'(W)SW)) W* -{n- {F'(W)SW) 

Qx(0,T) Jdnx(0,T) 

[ W* -(n- (F'(W)SW) = (44) 

This gives a method to evaluate 5 J but for each case the boundary conditions 



for W* are different because they must be chosen so that 5 J appears in (44) 



5.1 An example: the scramjet 

Consider a scramjet (see figure |4| with supersonic flow at the inlet and the 
outlet on which we observe the density at the outlet S and aim at tuning the 
parameters so that it be constant. 

When the flow is stationary and J is the integral on S X (0, T) of (p/poo — 
1) 2 /(2T) then, asymptotically in T, B = (1,0,0, ())//£, and 6 = 1. 

When the flow is hypersonic there is no condition on W at the outflow 
boundary S, therefore W a d is the whole space. Let W* on S be such that 

W* T n ■ {F'{W) = (BW - b)B on S (45) 
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With q 2 = u 2 + v 2 , S being parallel to the y-axis, this spells out as 



r 



w* 







(7-1) 2 



q — u 



1 ~~ \ 

(3-7)« (l-7)« (7-1) 

v u 



((7 - i)g 2 - JE) u (l + Ej - ( 7 -l)u 2 (1-7) 



vu yu 



= - - 1)(1 0) 

Poo Poo 

It is a non singular linear system for W* on S solvable as 



wt 



- £ -i) 

P 00 P 00 



7 -2_ 



7—2,-2 r 1 ■ 

u + v — jh u 



7-1 



-vWl 



w; = - 



-uWl 



7-1 

7-1 p 



/ 



(46) 



(47) 

(48) 
(49) 

(50) 



6 Numerical Flow Solvers 



In this section we wish to Implement numerically ( 47 1-( 50 ) , and compare with 



the conditions given by Automatic Differentiation of the Euler solver. 



6.1 Euler Flow Solvers 

The proposed method is a vertex-centered finite volume scheme applied to sim- 
plicial unstructured meshes which uses a particular edge-based formulation with 
upwind elements [lj [5] . The vertex-centered finite volume formulation consists 
in associating with each vertex of the mesh a finite volume cell, denoted 
Cj, which is built by joining each vertex to the middle point of its opposite 
edge. The common boundary between two neighboring cells C, and Cj is called 
dCij = dCi n dCj. An illustration of this construction in two dimensions is 
shown in Figure [3] 

The finite volume method assumes W piecewise constant on the cells and 
integrates ( [38] ) on each cell Cf. 

F(W l )-n t d 1 =\C l \ d ^+ V F\ ir f 11^7 = 0,(51) 
at Jar at ' p,r 

where is the outer normal to cell Cj, V(Pj) is the set of all neighboring vertices 
of Pi and F\ij represents the constant value of F(W) at interface dC^ from the 
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P;-side. The flow is calculated with a numerical flux function, denoted 

:= $y (Wi, W^ny) = F|y ■ / n, d 7 , (52) 

where iiy = j rij cfy. Several upwind numerical flux functions are available 



IdCi, 

and can be formally written: 

$y ( VF, , , ) = F ( w i) + F ( w i) . nii + d ( w. , , ny ) , (53) 

where the function d(Wi, Wj,n^j) contains the upwind terms and depends on 
the chosen scheme. Here the Roe approximate Riemann solver described in [53] 
is used. 



Roe's approximate Riemann solver is based on the Jacobian F'(W) : 

- ga±nS . niJ + |i (H ,, W y,^ (54) 

where A = P^T'F) evaluated for the Roe average variables W := (p, pu, pv, pE) 
as explained in 24 . For a diagonalizable matrix A = PAP -1 , |^4| stands for 
\A\ = P|A|P _1 . The eigenvalues of F'(W) are real and equal to u, u + c and 
u — c. 



Higher-order Finite Volume Methods. The previous formulation gives 
at best only a first-order scheme, higher-order extensions are possible with a 
MUSCL addition. The idea is to use extrapolated values Wij and Wji of W at 




Figure 3: Two finite volume cells Ci and Cj, and the upwind triangles Kij and 
Kji associated with edge PjP/. The common boundary dCij with the represen- 
tation of the solution extrapolated for the MUSCL correction is shown. 
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the interface dCij to evaluate the flux, cf. Figure [3j The following approxima- 
tion is performed: 

$ y = <!>,,( ir (/ . u;,,. n,, :• , 

with Wij and Wji which are linearly extrapolated by: 

W u = W k + ^ {VW) kl ■ Pk~P t with kl = ijorji, (55) 

where the approximate "slopes" (VW)ij and (VW)jj are defined on the edges 
by using a combination of centered, upwind and nodal gradients [T]. 

Such MUSCL schemes are not monotone. Therefore, limiting functions must 
be coupled with the previous high-order gradient evaluations to make sure that 



the scheme is TVD. Then VWm in (551 is substituted by a limited gradient 



denoted (VVF)'™ 1 . In this work we have used the three-entries limiter introduced 
by Dervieux which is a generalization of the Superbee limiter [1] . 



Boundary conditions. Slip boundary conditions are imposed for the walls 
because the flow is inviscid: U.n = 0. This boundary condition is imposed 
weakly by setting: 

$«*(Wi) = (0,j> i n i ,0) t . (56) 

For free-stream external flow conditions at inflow boundaries which ap- 
proximate infinity, a free-stream uniform flow Woo is known. There, the Roe 
flux is replaced by the Steger- Warming flux [2B] : 

$' x (Wi) = A + (Wi,n i )Wi + A-(Wi,n i )W 00 , (57) 

where A ± = ~(\A\ ± A) and A = F'(W). 

At free- stream outflow boundaries no condition is imposed so nothing is done 
to the finite volume scheme. 



Time discretization. Once the equations have been discretized in space, a 
set of ordinary differential equations in time is obtained: Wt — L(W) — 0. To dis- 
cretize in time, a high-order multi-step Runge-Kutta scheme is considered. Such 
time discretization methods, called SSP (Strong-Stability-Preserving), have non- 
linear stability properties which are particularly suitable for the integration of 
system of hyperbolic conservation laws where discontinuities appear. The opti- 
mal 2-stage order-2 SSP Runge-Kutta scheme introduced by Shu and Osher [25] 
is used in this study: 

W (1) = W n + StL(W n ), 

W n+1 = -W n + -W ( - 1) + -StL(W {1) ), 
2 2 2 x ' 



The scheme is stable with a CFL coefficient close to 1. Alternatively for (38 1 
an implicit time scheme can be used, such as Implicit Euler: 

\Y n+1 - W n 

|C<I St* +& i W n+1 = 0, (58) 
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where \tf (W n ) is the second order numerical flux associated with cell Ci given 
by the sum of the volume fluxes (53 1 and boundary fluxes. 



7 The Numerical Scheme for the Adjoint 

The adjoint of the discrete system is easier to establish when ((58|) is used. The 



linearized version of ( 58 ) is 



As the second order flux is too difficult to differentiate by hand, the second 
order Jacobian matrix is approximated by the first order one leading to: 



Thus the variation of the state variables 5W is given by 

A n SW n+1 = -^ 2 {W n ) . 



Consequently, following Section for Problem ( 38 1 , for a given functional 
J(W), the adjoint state W* is defined by 

A n *W* = J'(W) (59) 

where A n * is the transposed of matrix A" . Note however that it suffices to take 
A n * = {W n )) T because at convergence the solution does not depend on 
time and < $> 2 (W n ) tends to 0. 



7.1 Automatic differentiation of the numerical flux 

Tapenade |13j is a software which takes fortran or C programs as input and 
returns the differentiated programs; it is done by copy-paste of the program on 
the web page of the Tapenade internet site. One must specify the variable(s) 
with respect to which the differentiation is done. It works in direct or reverse 
mode and the later is very similar to the adjoint method. However the result may 
not be optimal. For instance if a program to compute the solution of F(u, a) — 
uses an iterative scheme like Newton's the differentiated code produced by 
Tapenade will also use Newton's method to compute u' a such that F' a (u : a) + 
F' u (u, a)u' a = even though this is a linear system which can be solved directly. 

In this study everything was differentiated by hand except the C-function 
which implements the numerical flux function which was passed to Tapenade 
for differentiation in direct mode. This way the full flux functions with MUSCL 
correction and slope limiter were differentiated. Alternatively we also wrote a 
solver where only a portion of the flux function is differentiated as explained 
below. Although full differentiation is safer, as we shall see, numerical test 
showed no difference when incomplete differentiation is used. 
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7.2 Manual approximate differentiation of the numerical 
fluxes 



One reason to prefer differentiation by hand of the flux functions of the finite 
volume method is complete control over the source code but also execution 
speed. 

Differentiation of Roe's flux In the Roe flux presented in Section [6] and 
given by (54) the Roe matrix (i.e., the dissipative matrix) is not differentiated: 



F , (W i )-im + \A(W i) W j 



A(Wi) + \A(W i ,W j )\) , 



and similarly for <&' w . Roe (Wi, Wj, ny). The analytic expression of is given 

in Appendix B. 



7.3 Differentiation of the boundary conditions. 

In order to differentiate the boundary conditions, we linearize the integrals which 
involve a boundary flux. Indeed, a Dirichlet condition does not require any 
computation as it fixes some state variable at a given value. 



The flux associated with the slip boundary condition is given by (561, hence 
using (76) in Appendix B, the following Jacobian is obtained: 



/ 



A 



Slip 



(7-1) 











\ 



\ 



q 

2 v 




o J 



where q 



2 i 2 

u + v . 





For the Steger and Warming type boundary conditions, where the flux is given 



by (57), the linearization is direct: 

A°°=A + {W i ,n i 



8 Numerical Results for a Scramjet with Out- 
flow Control 

A scramjet with supersonic flow at the inlet and the outlet has been computed 



with the flow solver described previously The discrete adjoint (59) was gener 



ated by hand using exact differentiated flux generated by Tapenade or manually 



using (60) 



When the flow is stationary and J is the integral on S x (0, T) of (p/poo — 
1) 2 /(2T) where S is the outflow boundary (parallel to the y axis), then we have 
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shown analytically that the boundary conditions at the outflow for the adjoint 
W* are given by ( 47 ) . 



8.1 Numerical Results 

The flow and the adjoint are shown on figure [4] 

The flow contains many shocks and some contact discontinuity. These em- 
anate from the rear tip of nozzle and the adjoint density is clearly discontinuous 
there. By contrast it is not discontinuous at the shock lines of the density. 

However shocks at the outflow boundary create discontinuous boundary con- 
ditions for the adjoints which, by analogy with Burgers' case, will propagate as 
discontinuities from right to left in the adjoint. The discontinuities are not on 
the flow shocks but around them on the characteristics which ends at the in- 
tersection of the shocks and the outflow boundary as in the one dimensional 
example of section [4] 



Boundary Conditions on the Outflow Boundary for the adjoint When 
the adjoint is generated automatically by A.D and hand transposition of the 
discrete equations we lose control over the boundary conditions at the outflow 
boundary, because are generated automatically and implement in weak form in 
the finite volume code. 

On ngure[5]the analytical values (47) are compared with the values generated 
by hand differentiation of the discrete solver. They agree quite well and there 
was no noticeable difference when the full Jacobian of the flux function computed 
by A.D. was used instead of the first order approximation. 

We noticed however that the shocks tended to be resolved better by the 
analytical expressions. Note also that W$ has a shock but no Dirac singularity 
in theory while the computed value seem to have a spurious Dirac singularity. 



9 Optimal Wing Profile with least Sonic Boom 

Consider an airfoil S at rest in a semi-infinite air domain f2. Air comes in from 
the left boundary L at supersonic speed and exits on the right by boundary 
R. Let the ground S be the lower horizontal part of the boundary ; thus the 
boundary of the domain is <9S! = SULUSUi?. The control criteria is on the 
pressure p: 

J=* J{p-Pq) 2 , therefore 8 J = J{p-p~o)5p. 
The aim is to find the best admissible £ to minimize J. 



9.1 Boundary Conditions 

Outflow Boundary Conditions on the Adjoint on R The flow is super- 
sonic so the situation is essentially the same as in section [8] for the scramjet 
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Figure 4: The scramjet: mesh (top), level lines of the density (middle) and of 
the adjoint density(bottom) 
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Figure 5: The scramjet, verification of (47 1. 
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except that J has no part on R so ( 46 ) has zero on the right hand side and both 
components of the normal n appear: 

W* ■F'(W).n = 

As F'(W).n is a non-singular matrix when the flow is hypersonic, this equation 
implies that all four components of W* are zero on R: 



W* = on R. 



(60) 



Inflow Boundary Conditions for the Adjoint on L The flow is supersonic 
on L so all components of the linearized flow 8W are equal to zero. All the 
characteristics are coming in so for the adjoint state the characteristics are 
coming out and hence 

W* needs no boundary condition on L. (61) 



Boundary Conditions for the Adjoint on the Ground S The shock is 
assumed to reflect perfectly so the vertical component of the velocity is zero: 
W 3 = 0. 

There is only one incoming characteristics for the linearized Euler system 
and for its adjoint as well. Assume that 6W3 — is the boundary condition 
for the linearized Euler system. In [8], Giles and Pierce gave a constructive 
argument to choose the boundary condition for the adjoint state, in this case 
W3 should be given. Let us verify directly that it is the right choice. 

On S the normal is (0, 1) T and the velocity is (it, 0) T , so 



n-F'(W) = 






(7-1), 






-(7 - l)u 



1 

(/. 




V 



rp 7 - 1 2 

u 

' 2 



7-1 
J 




Therefore 

W* ■ (n-F'(W))SW= (0 { -^u 2 5W! ~ (7 - l)uSW 2 + (7 - 1)5W 4 0)W* 

(62) 

On the other hand 



(7 - - u8W 2 + 5W 4 )W£ 



P = (7 - l )P e ~ 



7 - 1 (pu) 2 + {pvf _ W$ + W$ 



= (7-l)(W 4 - 



2Wi 



therefore 
By choosing 



2W? 



W? u 2 
5Wi - Trf L SW 2 + SWi) = (7 - 1)( T W - u6W 2 + SWi). 



W3 Is =P~Po, 
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(45) becomes (p — po)Sp and in the stationary regime (44) becomes 
SJ = 



dn\s 



W* ■ (n ■ F'(W)5W) 



(63) 



Boundary Conditions for the Adjoint on the airfoil E The flow is sub- 
sonic and tangent to E so u ■ n = 0. As for S one boundary condition is required 
for 8W . Let SW ■ n be given. An elementary calculus gives 

W*-(F'(W)-n)SW = {j-l)W* ■ n(^(u 2 + v 2 )SW! ~ uSW 2 - v5W 3 + 6W 4 ) 



-(w* + uw; + vw;)sw ■ n 



(64) 



As Giles and Pierce found [5] by another method, this expression contains no 
other <5VF-term but SW ■ n if and only if 



W* ■ n = on E. 
Proposition 1. Let W* be defined by 

d T W* + F'(W) T VW* = 0, W(T) = 0, 
W* ■ n = on E 

w*\ R = o, w;\ s =p~^ 

Then, asymptotically in time, 

8J=- I (W* + U ■ #2.3)^2.3 • n 



(65) 



(66) 



(67) 



where W2.3 is the vector (W2, Ws) T . 



Formula ( 67 ) is essential to set up descent algorithms for optimal shape de- 



sign to minimize J. The proof of the proposition is a consequence of ( 44 ) , ( 6 1 ) ( 63 ) , ( 65 ) . 
We have assumed that there was no shock on E; if it is not so, mean values have 
to be used for the integrand in (67). 

Remark 1. It may be more readable to rewrite (fi?| ) in terms of the components 
of the adjoint state (p* , (pU)* , (pe)*): 



SJ 



(p* +U ■ { P U)*)5{pU) ■ n 



(68) 



9.2 A Gradient Method for Shape Optimization 



Thanks to ( 68 ) we have now a method to modify the shape of the airfoil E to 
decrease J. 

Following [22j [20] , consider a variation E Q of E defined by the normal dis- 
placement a and given by 



E Q = {x + a(x)n(x) : x € E} 
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Denote W a the new flow variables and to simplify notations introduce V := pU 
By definition we have 

V(x) ■ n(x) = 0, V a (x + an(x)) ■ (x + an{x)) = 0, Vx G £ 

Assuming differentiability and ignoring higher order terms, the second equa- 
tion implies 

= V a {x) ■ n(x) + a ( (W(z) • n{x)) ■ n{x) + V{x) ■ n'(x) ) (69) 
where n' is the derivative with respect to a(x) of n^ a . The above is also 

dV 

SV ■ n\ s = V a (x) ■ n(x) = _a ("^ _ KV t) 

where t is the tangent vector associated with n - i.e. n rotated clockwise by 90 
degrees - because n' = —nt where k is the curvature (inverse of the radius of 
curvature) . 

Thus (see ( 68 1 ) by choosing 

a = -\(p* + u ■ (purx 9 ^- - K P u t ) 

for a small enough constant scalar A, J will decrease because, in this case 

5 J = -A Jjp* + U ■ (pUYfi 9 ^- ~ npU t f + o(X) 

This method was used by B. Mohammadi in [2U] and G. Roge et al in [6 
to reduce the sonic boom of a business jet by 30% without affecting the lift. 
Automatic differentiation was used, so we intend to reproduce the results with 
the code presented here in the future. Meanwhile, as for the scramjet, let us 
verify that the discrete adjoint converges to the continuous one. 



9.3 Numerical Tests 

A supersonic flow at Mach 2 is computed around a NACA0012 airfoil flying over 
ground at very low altitude, (see figure [6]). Euler flow and its discrete adjoint 
are computed on a flow-adapted mesh in 2D. The value of the adjoint on the 
ground S is compared with its analytical value. 

Results shown on figure [7] confirm that the discrete adjoint converges to the 
continuous one and that it is continuous at the shocks. 



Analysis of Figure [6] Here too we observe that the adjoint density is contin- 
uous at the shocks of p except near the ground. There the boundary condition 
on the adjoint forces a discontinuity which, by analogy with Burgers', should 
spread into shocks on the characteristics at the shock point. Since there is only 
one incoming characteristic we expect to have only one shock in the adjoint 
(while there was 2 in the case of Burgers'). On the reflected shock near the 
ground, the adjoint has a bump type variation. It is not clear whether it is 
physical or numerical. The adjoint also varies greatly near the NACA profile. 
A zoom shows that these variations are not discontinuities. 
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Figure 6: The NACA airfoil: the adapted mesh (top left),; the level lines of the 
density (top right) and (bottom, the figure is stretched horizontally) the level 
lines of the first component of the adjoint, i.e. the adjoint density. 
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'adj .dat' u 1 :4 + 
'adj.dat' u 1:6 x 



Figure 7: The NACA airfoil: comparison between , the component of the 
adjoint in duality with pv and p — po on S. Although the adjoint's equation is 
generated by automatic differentiation for the discrete case, it agrees with the 



continuous limit (see (66)) 



Discrete versus Continuous Adjoint (Figure [7]) For the minimizing of 
the sonic boom we have seen that the theory predicts that the boundary con- 
dition of the third component of the adjoint W 3 * is p — pq. On Figure [7] both 
quantities are displayed; the concordance is excellent. This is another indication 
that the adjoint obtained by automatic differentiation in reverse mode of the 
discrete solver converges to the continuous adjoint when the mesh is refined. 

Conclusion 

Industry is never happy to include an unreadable portion of code in their numer- 
ical tools; this is one reason to avoid using automatic differentiation to generate 
the adjoint of the Euler linearize system. By contrast A.D. is very practical and 
avoids code errors but theory doe not support yet the convergence when the 
mesh size decreases. 

On the whole this numerical study validates both the automatic generation 
of adjoints and the hand derivation with incomplete derivatives for the flux 
functions (i.e. without differentiation of the MUSCL and slope limiters). It 
is also reassuring to see on two examples that the discrete adjoints seem to 
converge to the continuous adjoints. 

One difficulty remains unanswered: whether the automatically generated 
adjoints are correct in the case of discontinuous boundary data. The theory 
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requires a precise boundary condition at the shock which the discrete adjoints 
don't seem to observe and even if they did a rather fine mesh would be needed 
to implement the condition anyway. However on the scramjet and on the airfoil 
this difficult did not show up in the numerical tests. 
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10 Appendix A: Calculus of Variation for Dis- 
crete Incompressible Euler Flows 

With reference to and Q, 

5J h = [ - u d )5u m+1 with 

n dt Jn 
= Y t J n Mx)l$<(y) - 6t8u%(x) ■ Vv%(y)]\ y = x - St<ix) dx (70) 

Now with a change of variable y = x — 5tu™(x) in the first half of the last 
integral, it is also 

1 

Jt 



v h (y + 5tv%{x))5v%{y)dy - / v h (x)8v%{x) ■ Vv%{x - 8tv%{x))dx (71) 
n+ Jn 



where fl + = {y := x — 5tu™{x) : x € f2}. 

The adjoint equation is obtained by replacing 6uh,8ph by Vh,qh and Vh,qh 
by u\,pt and adding the contribution of 5 J at time (m + l)5t on the right hand 
side: 

(-n+V<lh)u h - P h V-v h = (u h + -u d )v h 
n 01 Jn Jd 

+ Yt [ K m+1 (y + 8tu™ +1 (x))v h (y)dy - / v h u* h m+1 ■ V< +1 (y)dz (72) 
61 Jn+ Jn 

for all v h ,qh with v h \ r + = 0, v h ■ n\ r = 0. 

We choose to impose u^(T) — 0, u* h \ r + = 0, n£ • n|r = 0. 

Let us replace (w/ l7 g^) by (fc^, <5p?,,) where the tilde means that the value of 
Su™ +1 at the boundary nodes of T~ have been put to zero to be admissible: 

s j = Y,«t[ + v <w +1 K m - y prv ■ - / with 

' = £[/ ^ m+1 (x)^ +1 (j/)dy - ft / 5 M r +1 (^)«r +1 W ■ V^ +1 (y)dx] 



= y2l u h m (x)Su™(x - 8tu™{x))dx - ft / <5u™(s:)iift m (2;) ■ Vw™(y)da;] (73) 

It contains the equation of (Sun, Sph) with the test functions replaced by {u^,p* h ) 
it if it wasn't for the tildes. Let us call Sah the vector value P 2 finite element 
function which is equal to 8a at the nodes of and zero at other nodes. Since 
8uti = Suh + 8ah at all time steps, taking the boundary conditions into account 
we have up to 0(5t), o(5a,h) terms: 

SJ = StJ2 [ «™ ' « ■ V6a h ) + Sa h ■ (u* h m ■ Vv%)-p* h m V ■ 5a h ). (74) 
™ Jn 
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Appendix B: Jacobian matrix in 2D 

In two dimensions, the Euler equations could be symbolically rewritten: 

d t W + V • F{W) = , (75) 

where W = (p, pu, pv, pE) T is the conservative variables vector and the vector 
F represents the convective operator. The vector F may be decomposed as 

F(W) = Fi(W) e x + F 2 {W) e y with p = (7 - l)pE - l^l i^l + ^ an d 



F^W) 



( pu \ 
pu + p 

puv 

V ( P E+p)u J 



( 



and F 2 (W) 



pv 
puv 
pv 2 + p 



\ [pE+p)v J 

The Jacobian of p with respect to the conservative variable is: 



dp(W)^ (I 



-u , —V , 1 



where q 2 = u 2 + v 2 . After simplifications, the Jacobians read: 

1 



dF\jW) 
dW 



I 

(7-1) 2 2 



(76) 



-(7 - 3)u 



~uv 

2 



\ 

(7 — l)v 7 — 1 
u 



V ((7-1)9 ~lE)u 1 E- 



7-1 



(2u 2 + q 2 ) — (7 — l)uv 'yu J 



dF 2 (W) 
dW 





— uv 
(7-1) 2 2 

— q — v 



1 

u 



-(7 - l)u 



\ ((7 - l)g 2 - jE)v -(-y-l)uv ■yE - 1 -— - (2v 2 + q 2 ) jv ) 



-(7-3)w 



\ 


7-1 



For an edge e with n as normal vector, we have: ^(ly) :- 



( 

71 2 



dF(W).n 
dW 



q n x — u u • n u ■ n — (7 — 2)un x 
vn x — -yiun y 



^q 2 n y -vu- 11 



° \ 

71 n x 

JlUy 



u ■ n — (7 — 2)vn y 

y (71 q — 7-E-) u ' n (j~ + n x — 71 mu ■ n ^- + n y — -yivu ■ n 7U • n 
where u = (u,v), 71 = 7 — 1 and n = (n x ,n y ). 
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